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1.0  INTRODUCTION 


Ionograms  contain  as  basic  information  the  virtual 
range  versus  frequency.  The  more  advanced  digital  ionosondes 
also  measure  amplitude,  phase,  polarization,  incidence  angle 
and  spectral  signature  of  the  echoes.  This  allows  improved 
methods  for  the  conversion  of  the  ionogram  to  vertical  elec¬ 
tron  density  profiles.  This  report  discusses  two  different 
efforts  toward  that  goal.  Section  2  describes  the  automatic 
scaling  of  Digisonde  ionograms,  and  Section  3  investigates 
the  possibility  of  solving  the  Schroedinger  wave  equation  for 
different  potential  functions  in  an  effort  to  compare  the  WKP 
method  with  the  full  wave  solution. 


2.0  AUTOMATIC  PROCESSING  OF  DIGITAL  IONOGRAMS 


Automatic  processing  of  ionograms  has  become  possi¬ 
ble  with  the  availability  of  advanced  digital  sounders  [Bibl 
and  Reinisch,  1978],  Since  about  1970,  when  the  first  digital 
ionograms  were  routinely  recorded  on  magnetic  tape,  software 
and  hardware  techniques  for  the  automatic  extraction  of  the 
echo  traces  from  the  ionograms  [Bibl  et  al ,  1973]  have  been 
developed.  The  initial  concept  was  to  decide  immediately 
after  transmission  of  each  individual  frequency  which  of  the 
received  signals  are  echoes  and  to  retain  only  the  amplitudes 
and  heights  of  the  identified  echoes.  This  approach  is 
clearly  inadequate  for  many  scientific  investigations  in  a 
disturbed  ionosphere  where  the  researchers  want  to  see  the 
complete  range-versus-frequency  display,  i.e.  the  raw  ionogran 
with  all  the  nuances  in  the  signal  characteristics.  Bibl  an'’ 
Reinisch  [1978]  described  the  on-  or  off-line  printing  of  Digi- 
sonde  ionograms  using  an  optically  weighted  font  (Optifont)  to 
retain  the  digital  resolution  in  the  quasi-analog  ionogram 
display  (Figure  1),  establishing  the  transition  from  the  famil¬ 
iar  analog  to  the  digital  ionograms.  It  may  suffice  to  record 
on  magnetic  tape  a  limited  number  of  "identified"  echo  points 
per  frequency  for  the  archiving  of  routine  ionograms  or  for 
the  monitoring  of  ionospheric  trends  [Buchau  et  al,  1978], 

In  areas  with  frequent  occurrence  of  spread  F,  this  approach 
is  ineffective  unless  the  range  spread  for  each  frequency  is 
also  recorded. 

A  different  objective  is  the  automatic  calculation 
of  the  vertical  electron  density  profile.  In  that  case,  "the" 
vertical  echo  trace  must  be  extracted  and  oblique,  ducted  and 
multiple  echoes  disregarded;  a  trace  must  be  found  within  the 
spread  F  signals,  relying  on  amplitude,  Doppler  and  incidence 
angle  information  that  are  contained  in  the  digital  ionograms. 
This  is  only  possible  by  examining  the  ionogram  in  its  en¬ 
tirety.  This  report  describes  this  approach  for  the  auto- 
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Figure  1 


nat.ic  scaling  of  rigisonde  ionograms  recorded  at  the  Goose 
Bay  Ionospheric  Observatory  (53. 3N,  60.5V.7  geographic,  66. 6N, 

12. 1U  geomagnetic).  Once  the  ordinary  vertical-echo  trace  is 
extracted  from  the  ionogram,  the  electron  density  profile  can 
he  calculated  using  the  standard  lamination  technique  [Becher, 
1967J  or  the  profile-fitting  method  [Huang  and  Reinisch,  1982], 
the  latter  optimally  suited  for  the  case  of  automated  data. 

This  method  will  be  described  in  a  later  paper.  The  current 
report  discusses  the  accuracy  of  the  autoscaling  of  the  F- 
trace  during  disturbed  ionospheric  conditions,  e.g.  spread  F 
and  the  mid-latitude  trough.  Since  the  technique  for  the  auto¬ 
scaling  of  the  F-trace  is  not  quite  complete,  we  will  not 
describe  it  here. 


2.1  The  Ionogram  Data 

The  routine  ionograms  at  Goose  Bay  measure  and  re¬ 
cord  amplitude,  polarization,  incidence  angle  and  Doppler  fre¬ 
quency  as  a  function  of  frequency  and  range.  The  lower  part 
of  Figure  1  shows  a  quiet  amplitude  ionogram  containing  all 
signals.  Removal  of  the  non-vertical  and  y-polari zation  sig¬ 
nals  results  in  the  upper  ionogram  of  Figure  1,  which  is  much 
easier  to  scale  automatically.  As  will  be  shown  later,  the 
X-trace  data  are  not  discarded;  they  are  used  for  the  accurate 
determination  of  foF2.  For  bottomside  ionograms,  the  O-trace 
is  generally  better  presented  than  the  X-trace  and  our  auto¬ 
scaling  effort  concentrated  therefore  on  the  O-trace.  This  is 
in  contrast  to  topside  ionograms  where  0  and  X-polarization 
echoes  dominate  altematingly  [R.einisch  and  Huang,  19P2;  see 
also  the  Special  Issue  on  Topside  Sounding  and  the  Ionosphere 
of  the  Proceedings  of  the  IEEE,  Vol.  57,  No.  6,  June  19693. 

Ideal  ionograms  like  the  one  in  Figure  1  pose  no 
difficulties  for  automatic  scaling,  yet  they  are  useful  to 
illustrate  the  scaling  procedure.  To  start  the  scaling  of  the 
F-trace,  the  ionogram  for  h'  >  160  km  is  surveyed  and  the 


¥  — 


"center  window"  is  determined.  The  search  for  the  complete  re¬ 
trace  starts  at  the  center  window,  which  is  generally  the  best 
defined  part  of  the  trace  and  as  such  can  be  identified  most 
reliably. 

The  tape  recorded  Digisonde  ionograms  are  actually 
composite  ionograms,  since  each  pixel,  i.e.  frequency-range- 
bin,  is  accompanied  by  a  status  word  [Eibl  and  Reinisch,  1978] 
which  gives  the  signal  polarization  (C  or  X) ,  the  Poppler 
frequency  (four  spectral  lines),  and  the  approximate  incidence 
angles  (three  directions) .  The  current  report  discusses  the 
off-line  processing  of  ionograms  for  which  only  the  data 
stored  on  tape  are  available.  The  processing  algorithm  is 
tailored  toward  future  on-line  application  in  which  case  a 
larger  data  set  could  be  used,  like  separate  0  and  X  ionograms. 


2.2  Bottomside  Ionogram  Scaling  Algorithm-BISA 

Efficient  software  was  developed  that  determines 
the  O-echo  trace  h'(f)  in  about  five  seconds  of  CPU  time  on 
the  University  of  Lowell  Cyber  71  computer  or  the  AFC7  CPr  66 CO 
machine.  If  the  F-region  scaling  and  the  electron  density 
profile  inversion  routine  are  added,  a  total  of  about  ten  sec¬ 
onds  of  CPU  time  is  required. 

Figures  2a  and  2b  show  the  flow  chart  and  block,  dia¬ 
gram  for  the  main  BISA  program.  The  ionogram  data  are  entered 
from  magnetic  tape  and  unpacked  in  dual  arrays  of  128  ampli¬ 
tudes  and  128  status  words  for  each  ionogram  frequency.  Fach 
functional  unit  of  the  processing  sequence  is  briefly  described 
in  the  following  paragraphs. 

Noise  Elimination 

For  each  sounding  frequency  an  individual  noise 
threshold  is  determined  since  it  is  not  white  noise  but  man¬ 
made  interferers  that  generally  obscure  the  ionospheric  echoes. 
The  typical  signature  of  interference  is  a  fairly  constant 
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BISA  MAIN  PROGRAM 


Figure  2a.  BISA  Flow  Chart 
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distribution  of  amplitudes  over  the  range  bins,  while  the  echo 
amplitudes  generally  occupy  only  two  or  three  range  bins.  The 
most  probable  amplitude  found  over  all  range  bins  for  a  given 
frequency  is  used  to  determine  the  noise  threshold.  The  half- 
point  above  the  most  probable  amplitude  on  the  amplitude  prob¬ 
ability  curve  is  defined  as  signal  threshold.  Turing  severe 
spread  F  conditions  the  most  probable  amplitude  could  be  the 
signal  level,  setting  too  high  a  threshold.  To  avoid  this 
problen  two  probability  curves  are  formed,  one  for  the  lower 
64  and  one  for  the  higher  64  range  bins.  The  distribution 
that  peaks  at  a  lower  amplitude  sets  the  threshold  [Reinisch 
and  Smith,  1976]. 

0-Echo  Poles 

In  an  effort  to  minimize  the  required  CPU  time  the 
program  uses  for  certain  tasks  selected  data  points,  called 
the  poles.  For  each  frequency,  the  position  of  the  first  am¬ 
plitude  maximum  of  an  O-echo  between  two  zero  amplitudes  is 
defined  as  pole.  The  number  of  poles  per  frequency  is  not 
limited. 

Elimination  of  E  and  F  Fegion  Double  Echoes 

Before  beginning  the  trace  identification,  all  data 
points  that  are  double  echoes  are  eliminated.  An  F-region 
pole  at  twice  the  range  (±20  km)  of  a  lower  pole  is  checked  by 
comparing  the  average  amplitudes  over  nine  range  bins  of  the 
lower  and  upper  echoes.  If  the  upper  average  amplitude  is 
smaller  the  pole  and  all  poles  above  it  are  removed. 

The  same  method  is  applied  to  E-region  echoes,  ex¬ 
cept  that  a  5-bin  window  is  used.  Elimination  of  double  E 
or  Es  echoes  is  important  since  they  could  otherwise  be  mis¬ 
taken  for  F  echoes. 

F-Trace  Center 


During  disturbed  ionospheric  conditions,  it  is  often 
difficult  to  identify  the  main  overhead  echo  trace.  At  Goose 


Bay  this  is  especially  true  when  the  mid-latitude  F-region 
trough  moves  over  the  station.  The  direction  finding  capa¬ 
bility  of  the  Digisonde  identifies  some  of  the  oblique  echoes 
as  such,  but  not  all  of  them.  By  finding  the  frequency-range 
window  with  maximum  signal  energy  one  can  be  reasonably  sure 
that  the  trace  identification  starts  at  a  good  point  on  the 
main  echo  trace. 

For  each  range  the  pole  amplitudes  for  all  frequen¬ 
cies  from  the  first  ionogram  frequency  to  10  KIIz  are  added. 
Grouping  three  height  ranges  together,  a  search  is  made  for 
the  three  largest  peaks  in  this  sum  array.  A  5 -frequency  by 
3-range  bin  window  is  slid  along  each  of  the  three  maximum 
ranges,  and  the  sum  over  the  amplitudes  of  the  15  pixels  in 
each  window  position  is  formed.  The  window  position  with  tbe 
maximum  average  amplitude  is  designated  as  F-trace  center. 

Coarse  Baseline 

As  a  first  approximation  of  the  echo  trace  the  base 
line  is  constructed.  A  5-frequency  by  35-range  bin  window  is 
placed  over  the  trace  center,  and  the  median  range  h^  of  all 
poles  within  this  window  is  determined.  The  lowest  poles  in 
the  interval  h^  -  5  <_  h  +  20  bins  (bin  spacing  is  5  km) 

are  selected  as  baseline  points.  The  window  is  shifted  one 
step  toward  higher  frequency  and  the  next  baseline  point  is 
determined.  In  each  window  the  range  difference  D  between 
the  highest  and  lowest  pole  is  determined.  If  D  is  larger 
than  20  the  upper  height  of  the  window  is  increased  to  h  +  B 
If  no  pole  is  found  for  a  frequency  the  range  of  the  previous 
frequency  is  substituted.  The  process  stops  when  no  pole  is 
found  in  five  consecutive  frequencies.  It  is  assumed  that 
the  end  of  the  trace  is  found  and  the  baseline  point  is  giver 
the  range  128. 

Starting  again  .at  tbe  center,  the  process  is  re¬ 
peated  in  the  direction  of  lower  frequencies.  If  not  termin¬ 
ated  by  the  lack  of  proper  poles,  the  process  stops  at  f  = 
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foE  -  0.3  MHz.  Predicted  values  for  foE  are  tabulated  as 
function  of  time  of  day  and  season. 

Smoothing  of  Baseline 

The  sequence  of  heights  for  the  baseline  points  shows 
some  discontinuities  that  are  caused  by  high  interference  on 
some  frequencies,  and  blank  pixels  that  were  originally  occu¬ 
pied  by  oblique  or  X-polarization  echoes.  A  refined  linear 
smoothing  technique  is  applied  to  adjust  the  heights  of  the 
baseline  points. 

Consider  the  framed  trace  section  in  Figure  3.  The 
0-trace  points  are  marked  by  small  circles.  In  the  framed  sec¬ 
tion  the  original  baseline  points  are  shown  as  small  rectan¬ 
gles  if  they  differ  from  the  final  trace  points.  Smoothing  is 
performed  over  six  frequencies,  as  illustrated  in  Figure  A. 

The  shaded  areas  show  the  amplitudes  of  the  0-echoes,  the 
white  areas  the  amplitude  of  oblique  or  X-echoes .  The  height 
for  frequency  5.3  is  to  be  found;  the  heights  for  f  =  5. A,  5.5 
and  5.6  MHz  have  already  been  adjusted.  The  smoothing  starts 
at  the  trace  center  and  proceeds  to  the  lower  and  to  the  upper 
frequency  end,  respectively. 

A  straight  line  is  fitted  to  any  three  of  the  six 
frequencies,  resxilting  in  a  total  of  20  lines.  The  six  pixel 
amplitudes  for  each  line  are  added  and  the  line  with  the  maxi¬ 
mum  amplitude  determines  the  height  of  the  frequency  under  in¬ 
vestigation.  Only  three  of  the  20  possible  lines  are  shown 
in  Figure  A.  The  heavy  line  is  the  one  with  the  largest  am¬ 
plitude  sum  setting  the  height  for  f  =  5.3  MHz  to  360  km,  as 
marked  by  the  X. 

Determination  of  each  line  is  based  on  the  min-max 
approach.  The  equal-error  or  Chebyshev  [Scheid,  1968,  p.  269] 
line  is  determined  which  misses  three  given  data  points  by 
equal  amounts  with  alternating  signs  (see  the  deivations  -A, 

+A  and  -A  in  Figure  A).  There  exists  exactly  one  such  line 

h  =  c  +  df  (2.1) 
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TRACE  SMOOTHING  (LEFT  TO  THE  TRACE  CENTER) 


* 


which  is  defined  by  the  following  three  equations: 

c  +  df^  -  h^  *=  A 

c  +  df2  -  h2  “  -A 

c  +  dfj  -  h-j  *  A 


The  solutions 
c  ■ 

d  - 


for  c,  d  and  A  are: 


h^  +  A  -  dfj 


h3"h 


1 

1 


A  = 


)]. 


(2.D 

(2-3) 

C.'O 


C  ) 

r  » 

C  .  ) 


Hyperbolic  Trace  Fitting  in  Cusp  Region 

If  the  ionogram  data  have  polarization  tagging,  sorr.r 
curve  fitting  can  be  performed  in  the  h'-f  domain.  The  cusp 
region  of  the  ionogram  trace  can  be  described  with  reasonable 
accuracy  by  a  polynomial  or  a  rational  function.  A  simple 
rational  function  like 


K 


1 

T 


b  <  0 


(2.8) 


has  a  singularity  at  f  *  -  a/b  which  serves  well  to  fit  the 
trace  near  the  critical  frequency.  The  same  function  trans¬ 
posed  to  lower  frequencies  by  half  the  gyrofrequency  fH  de¬ 
scribes  the  corresponding  O-trace: 

ho  =  h0  +  F+^bf*  ;  f*  =  f  +  7  fH300  (2.9) 


A  reference  height  of  300  km  is  used  for  the  calculation  of 
fH.  Fitting  simultaneously  two  curves  to  the  cusp  data  makes 
the  method  very  robust.  Even  if  one  of  the  traces  is  absent, 
the  critical  frequency  is  determined  from  the  other  trace. 

The  hyperbolas  with  the  maximum  sum  of  amplitudes  in  narrov7 
strips  around  them  are  selected  to  represent  the  0  and  X 


< 
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(  t  the  critical  frequencies.  The  optimum  positioning 

of  the  hyperbolas  is  illustrated  in  Figure  5  for  a  Goose  Bay 
ionogram  with  spread  F.  The  range  interval  over  which  the 
amplitudes  are  summed  is  given  for  each  frequency  hy  the 
ranges  of  the  preceding  and  the  following  frequency  as  shewn 
in  Figure  5.  The  Digisonde  ionogram  used  in  this  figure  does 
not  show  the  signal  amplitudes  but  instead  the  polarization 
tags  0  and  X,  and  the  letter  E  indicating  oblique  echoes. 

Only  O-amplitudes  contribute  to  the  O-hyperbola 1 s  amplitude 
summation  and  X-amplitudes  to  the  X-hyperbola.  The  hyperbolic 
traces  for  the  cusp  region  that  the  program  determined  for  the 
ionogram  in  Figure  5  are  marked  by  heavy  dots  (O-trace)  and 
X's  (X-trace).  The  hyperbolic  section  of  the  trace  ends  where 
the  X-hyperbola  intersects  the  baseline  0.5  MHz),  and  the 
baseline  then  defines  the  trace  toward  lower  frequencies. 

Equations  2.8  and  2.9  have  three  coefficients  hp,  a, 
and  b,  which  are  determined  according  to  the  method  of  min-nax 
rational  functions  [Scheid,  1968,  p.  289],  The  baseline  is 
used  as  a  guide  in  finding  the  best  fit.  Starting  at  the  high 
frequency  end  of  the  baseline,  a  hyperbola  is  fitted  to  the 
last  three  baseline  points  by  minimizing  the  distance  x  from 
the  hyperbola  (see  Figure  6).  The  hyperbola  is  determined  in 
such  a  way  that  the  distances  are  x,  -x  and  x,  where  x  is  a 
positive  or  negative  number.  The  range  deviations  of  the 
three  points  (f^,  hj)  from  the  curve  are 

hi  -  {h0  +  rTT¥T}  "  x  i  -  1,  2,  3.  (2.10) 


These  equations  can  be  written  in  the  form  (deleting  the 
primes  on  the  heights): 


(hi 

-  hp  -  x)a  + 

(h^  -  hQ  -  x) f ^ 

O 

II 

rA 

1 

(2.11) 

(h2 

-  hg  +  x)a  + 

(h2  -  hc  +  x)f2 

b  -  1  «  0 

(2.12) 

(h3 

-  tiQ  -  x)a  + 

(h3  -  hQ  -  x)f3 

b  -  1  -  0 

(2.13) 
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o  +  bf 


From  (2.11)  and  (2.12)  one  obtains 


(2.14) 

(2.15) 


Elimination  of  a  and  b  leads  to  a  quadratic  equation  for  the 
determination  of  x  CScheid,  1968,  p.  289]: 

ih1  -  hQ  -  x  O'!  -  h0  -  *>  fi  1 

jh2  -  hp  +  x  (h2  -  hp  +  x)  f-2  1  =  0 

I 

i  h3  •  Hq  *  x  (h^  "  hp  •  x)  f  3  1  (2 . 16) 

or: 


Ax2+Bx+C=0 
The  coefficients  A,  B  and  C  are: 

A  =  ?  (fl  -  f3) 

B  =  Bq  +  AhQ 
c  =  CQ  -  ci  hQ 

where 

B0  =  -  f1(h3+2h1-h2)  -  f 2 (h3~h1)  -  f3(h2-2h3-h1) 
C q  *  f  l.hi (h3“h2)  "b  f ph2  (h^-h3)  +  f3h3(h2~h^) 

Cl  *  f 1 (h3”h2)  +  f2 ^h^~h3)  +  f3 Ch2-h3) 


(2.17) 

(2.18) 

(2.19) 

(2.20) 

(2.21) 

(2.2?) 

(2.23) 


Selecting  the  smaller  of  the  two  solutions  for  equation  2.17 
we  get : 


The  coefficients  E  and  C  are  functions  of  the  base  height 
(equation  ?..8),  which  is  not  known.  Inspection  of  equation 
shows  that  \r.  j  will  be  a  minimum  if  4AC/E?-  0.  This 
nay  not  ho  possible  to  achieve',  hut  one  can  fry  to  find  the 
minimum  value  of  |AC/B2  |  .  This  is  done  by  setting  the  deriv¬ 
ative  with  respect  to  hp  equal  to  zero: 


3FT  <->  “  0 


0  B2 
resulting  in 


(2.25) 


h0  = 


CL  B0  +  2  AC0 

- IT, - 


(2.2.6) 


Once  hg  and  x  are  calculated  from  (2.26)  and  (2.24)  the  coeffi¬ 
cients  a  and  b  can  be  calculated  specifying  the  X-hyperbcla. 

The  corresponding  O-hyperhola  is  simply  determined  by  replac¬ 
ing  f  by  f  .  Finally  the  0  and  X  sums  are  determined. 

It  is  of  course  unlikely  that  the  hyperbola  fitted 
tc  the  last  three  points  (highest  frequency)  of  the  baseline 
will  be  a  good  approximation  for  the  X- trace  in  the  cusp  re¬ 
region.  So  the  process  is  repeated  for  the  next  three  points 
and  so  on,  dovm  to  frequencies  at  the  center  window,  and  then 
again  for  non-abutting  points.  If  v.Te  number  the  baseline  data 
points  from  1  to  K,  we  can  write  the  triple  point  sets  to 
which  hyperbolas  are  fitted  in  the  following  way: 

GROUP  1  GROUP  2  GROUP  3 


rl 

P2 


rl 

Po 


PN-2  PK-1  PN  PN-4  PN-2  PN  PN-6  PU-3  PH 

The  total  number  of  trial  hyperbolas  is  close  to  2N  requiring 
precious  CPU  time.  A  factor  of  two  in  time  could  be  gained  by 


deleting  Group  1  of  the  guiding  points.  From  all  these  hyper¬ 
bolas,  the  one  with  the  maximum  0  and  X  amplitude  sum  is  se¬ 
lected  for  the  trace  representation. 

0-Echo  Trace 

Combining  the  O-hyperbola  with  the  smoothed  baseline 
gives  the  echo  trace  for  the  O-polarization  echoes.  If  the 
Digisonde  outputs  complete  0  and  X  ionograms,  as  can  be  done 
with  the  new  Digisonde  256,  one  can  devise  a  complete  C  and  X 
trace.  This  would  be  of  interest  mainly  for  stations  whe  want 
to  study  the  ionization  valley  between  E  and  F-region. 

Figures  7  and  8  illustrate  the  performance  of  the 
scaling  algorithm  for  Goose  Bay  ionograms  with  light  spread. 
The  preliminary  results  for  the  E-region  scaling  are  included 
in  these  figures,  together  with  the  corresponding  electron 
density  profiles. 

Parameter  Extraction 


At  the  present  time,  the  program  determines  the  fol¬ 
lowing  F-region  parameters:  critical  frequency  of  the  F- 
region,  f oF2 ,  minimum  height  h'F,  the  11(3000)  factor  and  the 
spread  F  figure.  The  foF2  value  is  obtained  in  the  hyperbolic 
fitting  procedure  as  -a/b.  The  smallest  range  value  in  the 
F-trace  is  used  as  h'F.  The  M(3000)  factor  is  obtained  from 
the  transformed  oblique  ionogram  which  is  calculated  from  the 
vertical  ionogram  trace  by  multiplying  each  frequency  with  the 
transmission  factor  M(h')  [Smith  et  al,  1979,  p.  123: 


£ob  =  M(h’).fv 


(2.27) 


The  transmission  curve  M(h')  is  calculated  by  fitting  a  poly¬ 
nomial  to  the  UPSI  specified  data  set  [UP.SI  Handbook  of  Iono¬ 
gram  Interpretation  and  Reduction,  Second  Edition,  November 
1972,  p.  21;  World  Data  Center  A  Report  UAG-233.  The  MUF(300G), 
obtained  as  the  highest  frequency  in  the  oblique  ionogram,  is 
divided  by  foF2  to  find  the  hf(3000)  factor. 
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determined 


For  each  ionogram  an  average  spread  F  figure  SF 
by  calculating  the  mean  range  spread: 


is 


SF 


1 

ToTT^TmTnT 


loF? 


I 

fminF 


Ah(f)  . 


(2.28) 


The  upper  limit  of  the  range  spread  Ah(f)  is  defined  by  a  drop 
in  average  echo  amplitude  (averaged  over  four  range  bins  or 
20  km)  by  12  dB  or  more  below  the  average  trace  amplitude. 


2.3  Performance  Evaluation  of  the  Ionogram  Scaling  Algorithm 

It  is  always  easy  to  automatically  scale  a  few  se¬ 
lected  ionograms  by  tailoring  the  scaling  algorithms  to  these 
specific  ionogram.  The  real  test  is  the  application  of  the 
algorithm  to  a  large  number  of  ionograms  covering  day,  night 
and  twilight  as  well  as  quiet  and  disturbed  conditions.  To 
test  the  performance  of  the  PISA  program  one  month  of  iono¬ 
grams,  i.e.  approximately  2000  ionograms,  for  January  1980 
from  Goose  Bay,  Labrador,  were  processed  and  the  parameters 
of  the  approximately  700  hourly  ionograms  were  compared  with 
the  manually  scaled  values. 


2.3.1  Results  of  the  Autoscaling 

Figures  9,  10  and  11  show  the  diurnal  and  day  to  day 
variations  of  foF2,  h'F,  and  the  spread  F  figure  for  January 
1980  using  the  Optifont  to  visually  enhance  the  variations. 
With  three  ionograms  per  hour  the  72  daily  values  are  printed 
in  one  line  with  consecutive  days  following  each  other.  The 
foF2  values  in  Figure  9  are  printed  in  increments  of  0.5  MHz 
although  the  original  scaling  resolution  is  0.1  MHz.  From  a 
pre-sunrise  minimum  below  A  MHz  at  07  AST,  foF2  increases 
rapidly  after  sunrise.  When  foF2  increases  from  11.5  to  12 
MHz  the  Optifont  numbers  go  from  15  to  0,  so  the  12  MHz  con¬ 
tour  line  in  Figure  9  is  indicated  as  the  sharp  break  in  in- 
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tensity.  The  highest  foF2  values  were  observed  around  January 
when  foF2  reached  values  between  14.5  and  15  MHz.  Later  in 
the  month  foF2  stayed  at  or  below  12  MHz.  The  occurrence  of 
the  F-layer  trough  becomes  visible  as  a  sudden  decrease  in 
foF2  during  the  early  evening  hours  on  days  3,  5,  7,  11,  127, 

28  and  29.  Days  12  and  13  show  irregular  foF2  values  indicat¬ 
ing  ionosphericallv  disturbed  conditions.  Indeed,  the  magnetic 
A-index  was  above  45  on  these  days,  while  A  <  25  on  the  other 
days . 

Figure  1C  displays  the  h'F  values  in  10  km  incre¬ 
ments  using  the  same  format  as  the  previous  figure.  The 
height  variations  during  the  day  hours  are  small  except  during 
the  disturbed  days.  The  disturbed  days  and  the  occurrence  of 
the  trough  are  even  better  indicated  by  tbe  spread  F  figures 
displayed  in  Figure  11.  The  numbers  printed  give  the  average 
range  spread  in  multiples  of  5  km. 

The  difficulty  of  the  scaling  task  is  illustrated 
by  the  sequence  of  ionograms  in  Figure  12  recorded  in  Coose 
Bay  on  7  January  1980  between  1819  and  2039  AST.  It  is  the 
time  when  the  F-region  trough  moves  over  the  station  causing 
spread  F  and  oblique  echoes.  The  results  of  the  autoscaling, 
indicated  by  the  small  filled  circles,  are  in  good  agreement 
with  manual  interpretation.  The  amplitude  ionograms  in  Figure 
12  contain  the  0  and  X-echoes  as  well  as  oblique  returns.  If 
reflections  from  the  trough  walls  or  other  irregularities, 
libe  at  1959  AST,  produce  several  retardation  cusps  in  the 
ionogram,  the  program  most  libely  selects  the  one  with  tbe 
strongest  amplitudes.  A  manual  evaluator  may  scale  foF2  at 
1959  AST  as  2.7  or  3.2  MHz,  while  the  program  scaled  4.3  MHz 
and  jumped  to  the  lower  value  at  2019  AST.  These  differences 
in  "interpretation"  are  the  main  contributor  to  errors  in  foF2 
above  0.5  MHz. 


2.3.2  Comparing  Manual  and  Autoscaling 


Since  manually  scaled  parameters  for  the  hourly 
ionograms  were  available,  they  were  used  as  reference  against 
which  the  autoscaled  values  are  compared.  For  a  total  of 
about  570  ionograms  the  parameters  foF2,  h'F  and  fminF  were 
investigated  for  their  accuracy.  Ionograms  with  technical 
errors  were  removed  from  the  data  base. 

Figure  13  shows  the  error  distribution  function  of 
foF2  where  the  error  is  defined  as  the  difference  between  the 
manual  and  autoscaled  value.  The  error  distribution  is  Gaus¬ 
sian  in  shape;  82%  of  all  ionograms  have  errors  of  less  than 
C.5  MHz,  and  93%,  fall  within  the  1  MHz  limit.  If  the  analysis 
is  limited  to  ionograms  without  spread  F  the  statistics  are 
only  slightly  improved  as  can  be  seen  in  Figure  14.  This 
demonstrates  EISA's  success  in  the  scaling  of  disturbed  iono¬ 
grams.  The  diurnal  variation  in  the  autoscaling  of  foF2  can 
be  checked  in  Figure  15.  It  is  apparent  that  the  scaling  accu¬ 
racy  is  higher  for  the  day  ionograms  than  for  the  night  mea¬ 
surements.  This  is  expected  because  of  the  complicated  signa¬ 
ture  of  the  night  ionograms  at  Goose  Bay. 

Comparison  of  the  median  foF2  values  shows  very  good 
agreement  between  manual  and  autoscaling  (Figure  16).  Two 
larger  deviations  occur  at  1700  and  .1900  AST  which  are  caused 
by  differences  in  interpreting  the  ionograms  in  the  presence 
cf  oblique  signals.  The  manual  median  values  are  based  on  the 
URSI  conventions  regarding  qualifying  letters,  no  qualifica¬ 
tions  were  used  for  the  autovalues.  This  explains  some  of  the 
smaller  deviations. 

The  error  distribution  functions  for  fminF  and  h'F 
(Figures  17  and  18)  verify  the  excellent  performance  of  BIS/. 
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(USING  580  IONOGRAMS,  JANUARY  1980,  GOOSE  BAY,  LABRADOR) 


2.4  Summary  and  Outlook 

By  comparing  the  autoscaled  ionospheric  parameters 
with  manually  scaled  values  for  a  whole  winter  month  at  Goose 
Bay,  it  could  be  shown  that  the  BISA  program  can  be  success¬ 
fully  employed  even  at  stations  with  an  unusual  high  rate  of 
ionospheric  disturbances.  The  program  must  now  be  completed 
by  adding  tbe  E-region  scaling  algorithm  and  by  specifying 
some  qualifying  parameters.  The  basic  concepts  for  these 
tasks  have  been  developed,  so  that  the  program  can  be  com¬ 
pleted  within  the  next  year.  Solutions  for  the  complementary 
task  of  converting  the  autotrace  into  an  electron  density  pro¬ 
file  are  based  on  the  profile-fitting  method  [Kuang  and 
Reinisch,  1982],  Ve  originally  developed  and  applied  this 
method  to  autoscaled  topside  ionograms.  The  necessary  adjust¬ 
ments  for  bottomside  profiles  will  be  made  in  the  next  year. 

Efficiency  in  the  computer  coding,  one  of  the 
achieved  design  goals,  makes  it  possible  to  implement  the 
software  in  an  on-line  ionogram  processor,  the  so-called  real¬ 
time  ionogram  scaler  (RIS) .  The  RIS  will  be  based  on  the  8086 
microprocessor  and  uses  Fortran  86  as  source  language. 


3.0  COMPARISON  OF  FULL  WAVE  THEORY  AND  WKB  APPROXIMATION 
FOR  DETERMINING  IONOSPHERIC  DENSITY  PROFILES  FROM 
REFLECTION  COEFFICIENTS 


3.1  The  Full  Wave  Theory  and  the  WKB  Approximation.  Time- 

Delav  From  the  Phase  of  the  Reflection  Coefficient. 

[Cohen  et  al,  1982] 

3.1.1  The  Basic  Wave  Equation 

Let  E(z,t)  be  a  component  of  the  electric  field  per¬ 
pendicular  to  the  z-axis  or  vertical  axis  along  which  the 
electron  density  is  stratified.  Then  with 

Ek(z,t)  =  e"iwt  Ek(z)  , 

(iu  =  ck)  (3.1) 

the  equation  for  Ej.(z,t)  is  [Budden,  I960,  page  129] 

E,(z)  +  k2  n2  E,(z)  =  0  .  (3.2) 

dz2  k  k 

In  (3.1)  u  is  the  angular  frequency,  k  =  ^  is  the  wave  number 
in  free  space  and  n  is  the  index  of  refraction  given  by 

n 2  =  1  —  ^-(-2)ZfTre2  (3.3) 

mc2k2 

where  m,  e,  c  are  the  mass  of  the  electron,  the  charge  of  the 
electron  and  velocity  of  light  in  Gaussian  cgs  units. 

In  our  discussion  we  have  ignored  electron  colli¬ 
sions  and  the  effect  of  the  earth's  magnetic  field.  This 
physical  (as  opposed  to  mathematical)  approximation  is  often 
assumed  in  ionospheric  sounding  methods.  At  the  magnetic 
equator,  the  earth's  magnetic  field  is  parallel  to  the  earth's 
surface  and  hence  the  polarization  of  the  transmitted  wave  can 
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be  chosen  so  that  the  magnetic  field  does  not  affect  the  scat¬ 
tered  wave.  In  other  situations,  the  effect  of  the  magnetic 
field  may  also  be  eliminated  [Budden,  1966].  The  effect  of 
collisions  is  negligible  at  higher  altitudes  (z  >  100  km) .  One 
of  the  ultimate  objectives  of  the  present  train  of  research  is 
to  determine  whether  these  physical  approximations  lead  to 
errors  greater  than  those  made  by  the  mathematical  approxima¬ 
tions.  In  any  case,  we  shall  assume  as  is  customary  that  the 
use  of  Eq.  (3.3)  for  the  index  of  refraction  is  a  good  approxi¬ 
mation  for  use  in  ionospheric  sounding  in  certain  geographical 
regions . 

Equation  (3.2)  can  also  be  written 

—  E.  (z)  +  (k2  -  V (z))  E,(z)  =0  (3. A) 

dz2  k  k 

where  V(z)  is  given  by 

V (z)  =  K  N(z),  with  K  =  .  (3.Aa) 

me2 

Equation  (3.4)  is  the  one-dimensional  Schroedinger 
equation  which  has  been  exhaustively  studied.  The  potential 
V(z)  is  essentially  the  number  density  N(z). 

In  the  direct  problem  of  reflecting  electromagnetic 
waves  from  the  ionosphere  we  assume  V(z)  (or  N(z))  is  very 
small  for  z  <  Zq  and  z  >  z^. 

We  look  for  solutions  of  (3. A)  which  behave  like 
E^(z)  =  eiks  +  b (k)  e"ikz  for  z  <  zQ 

=  t (k)  eikz  for  z  >  z1  .  (3.5) 

The  quantities  b(k)  and  t(k)  are  called  the  reflec¬ 
tion  and  transmission  coefficient  respectively. 
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For  a  wave  with  wave  number  k ,  we  have 


rcz.t)  -  e1H*-ct>  +  b(k)  e-lt<:!!+ct),  z  <  z0 

-  too  eik(2"ct),  z  >  Zj  .  (3.0.) 

The  boundary  conditions  (3.6)  are  interpreted  to 

2  if 

mean  that  a  plane  wave  with  wave  number  k  =  y—  moves  initially 
toward  the  scattering  potential  V(z)  and  then  is  partially  re¬ 
flected.  A  transmitted  wave  on  the  other  side  of  the  poten¬ 
tial  propagates  in  the  Z-direction. 

Actually,  one  never  sends  an  infinite  plane  wave  to¬ 
ward  the  potential  (or  ionosphere).  Instead  one  sends  a  pulse 
containing  several  wave  lengths.  This  pulse  can  be  represented 
by 

E(z.t)  -  /  A (k)  e-ia;t  Fk(z)  dk  ,  (3.7) 

—  00 


and  is  thus  a  superposition  of  the  plane  waves.  The.  amplitude 
factor  A(k)  has  its  peak  value  near  or  at  the  value  of  k  =  k^, 
where  kp  is  the  mean  wave  number  which  appears  in  the  pulse. 
One  sends  a  pulse 


^incident  ^ 


/  A(k)  eik(2-ct)  dk 


O.f) 


and  gets  back  a  reflected  pulse 

FrefUc«d<2't>  *  /"A(k)  b(k>  e'lk<2+ct>  dk.  <3.") 

a*  on 

Since  Etocldent(*.t)  and  *reflected<z. O  are  of  finite  extent 
we  can  ask  for  the  time  for  the  reflected  wave  to  return  to 
the  transmitter.  This  time  will  depend  on  kp.  It  is  given 
by 

T(k0)  -  £  evaluated  at  k  -  k0  (3.10) 
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where  <?(k)  is  the  phase  of  b(k);  i.e. 

b(k)  -  |b(k)i  e1*00  .  (3.11) 

The  virtual  height  h  is  given  by 

hv  =  i  cT(k0)  .  (3.11a) 

It  should  be  remembered  that  the  theory  thus  far  is  exact. 

(In  using  (3.10)  for  the  time  delay  it  is  convenient  to  thinV 
of  the  transmitter  as  being  located  at  z  =  0.) 


3.1.2  The  Approximate  Inverse  Problem 

The  direct  problem  is,  for  our  purposes,  the  follow¬ 
ing:  Given  V(z)  (or  equivalently  N(z)),  find  T(k). 

The  inverse  problem  is:  Given  T(k)  for  all  k,  find 
V (z)  (or  N(z)) . 

The  approximate  solution  is  obtained  from  the  1YP 
approximation.  Assume  Ek(z)  has  the  form 

Ek(z)  *  A  eiF(z)  .  (3.12) 

Substitute  (3.12)  into  equation  (3.2)  to  obtain 

(d|)2  =  k.2  n2  +  i  (3.13) 

which  becomes  a  Riccati  equation  if  y  =  dF/dz.  An  iteration 

procedure  can  be  started  by  assuming  on  the  right  of  Fq.  (3.13) 

d  2  F 

that  k2n2  is  large  compared  to  ^r»  i.e.  n/A  is  large  or  n 
varies  slowly  compared  to  the  wave  length  [see  Kay,  1954  for  a 
more  careful  analysis  of  domain  of  validity].  The  first  iter¬ 
ation  gives  F'  «  ±  kr  while  the  second  yields 

F*  -  ±  kn  (1  +  i£l-)1/2  :  ±  kn  (1  +  ----■)  *  ±  kn  + 

k2n2  2k2n2  "n 
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which  gives  F  by  a  integration.  Since  Eq.  (3.13)  is  non-linear, 
the  two  particular  solutions  of  (3.13)  cannot  be  superimposed 
to  find  a  general  solution,  however,  the  solution  F  can  be 
substituted  into  Fq.  (3.12)  to  obtain  particular  solutions 

z 

Ek(z)  =  A  [n(z)]“1/2  exp  [*  i  k  /  n(z')  dz']  (3 . 1 A> 

z0 

to  Eq.  (3.2).  It  should  be  remembered  that  n(z)  is  a  function 
of  k  (see  Eq.  (3.3)).  We  have  suppressed  this  dependence  on 
the  notation  for  simplicity.  If  one  particular  solution  to 
Eq.  (3.2)  is  known,  it  can  be  used  to  give  a  first  order  lin¬ 
ear  equation  which  can  be  solved  by  quadratures  thereby  giving 
the  general  solution.  Alternatively,  two  distinct  particular 
solutions  can  be  added  to  give  the  general  solution  of  Eq. 

(3.2)  which  is  (within  this  approximation) 

Ek(z)  =  A(k)  Cn(z)]"1/2  exp  [+  i  k  /  n(z’)  dz1]. 


2 

+  B(k)  [n(z)]"1/2  exp  [-  i  k  /  n(z')  dz'].  (3.15) 

z0 

Using  the  first  boundary  condition  of  Eq.  (3.6)  end  noting 
n(z)  =  1  for  z  «  Zq,  one  finds 


Ek(z)  =  eikz  +  b(k)  e~ikz  (z  <  zQ) 

-  A(k)  e'ik2°  elkz  +  E(k)  e^0  . 


Hence  A(k) 


ikz 

e 


0 


B(k) 


ikz 

b  (k)  e 


(3.16) 


Let  zk  be  the  value  of  z  for  which  n(z)  *  0,  i.e.  from  Eq. 
(3.3)  and  (3.4a) 

V(zk)  -  k2  .  (3.17) 
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Our  picture,  the  usual  one,  is  that  V(z)  =  0  for  z  <  Zq  and 
increases  monotonically  toward  a  maximum  as  z  increases  above 
Zq.  As  k2  varies  from  0  to  the  maximum  of  V(z^,),  a  z^  is  de¬ 
fined  by  (3.17). 

From  (3.15),  E^(z^)  -*■  »  unless 

zk  zk 

A(k)  exp[i  k  /  n(z')dz']  +  E(k)  exp[-  i  k  /  n(z')  dz ' ]  =  0. 


Thus 


B(k)  -  -  A(k)  exp[2ik  /  n(z')  dz '  ] 


(3.18) 


ikz  x 

-  e  ®  exp[2ik  /  n(z')  dz '  ]  . 

zrt 


(3.19) 


Or  from  (3.16) 


2  ^ 

b(k)  =  -  e  “  0  exp [ 2 ik  /'  n(z')  dz']  . 


(3.20) 


The  phase  <f> (k)  of  Eq.  (3.11)  is 


(k)  -  7r  +  2  k  [zn  +  /  n(z')  dz']  . 


(3.21) 


zk 

2k  r  8 


T (k)  “  f  Czo  +  /  n<2’)  dz’^  +  —  /  ?£  n(z')  dz' 


+  rL-><V  irr  • 


(3.2?) 


But  n(z,)  -  0  by  definition  (see  Fqs.  (3.3)  end  (3.4)),  also 


Sn(z)  8  1  JL  2  t ,/  \ 

-JT1  ■  57  IT  A  -v<z> 


i _  J\r  2  If  f  A.  _ 1 


-  —  /k2-V(z)  + 
k2 


✓k2-V (z) 


*  k  n(z)  +  IcnTzT  * 


Thus  finally 


(3.2?) 


2  r_  -L  fk  dz’ 


T<k>  -  I  Cz0  +  l  nTFT-1 


TOO  =  |  Czp  +  k  / 


(3.24) 


This  is  a  well-known  integral  equation  for  V(z).  The 

c.Z  r\ 

„  4-1 _ _  J  ^1.4.  w  _ _ _ I _ 4.1*  ^  » _ .T _ J.V  _  _  *  _ _ 1 


first  term  on  the  right 


represents  the  time  for  the  signal 


to  reach  Zq  and  return  to  the  transmitter  at  z  *=  0.  If  k^  is 
the  lowest  value  of  k  for  which  there  is  a  reflection,  then 
approximately 


cT^) 


(3.25) 


Now  that  Zq  is  determined  approximately,  one  can  find  V(z) 

(or  N(z))  using  Abel's  method  of  inversion  [Faddeev,  1958], 

For  each  value  of  k  one  gets  T(k)  experimentally.  One  then  can 
find  V(z)  (Zq  <  z  <  z^.)  using  Abel’s  method.  This  is  the  tisual 
procedure.  Note  if  k2  >  V  (z) ,  the  method  fails. 


3.1.3  The  Exact  Inverse  Problem 

The  exact  or  full-wave  theory  version  of  the  inverse 
problem  is  developed  in  Kay,  1954  and  Kay  and  Moses,  1956. 


Let  b(k)  be  the  reflection  coefficient  as  in  Eq. 

(3.6).  The  potential  can  be  recovered  from  a  knowledge  of 
b(k)  through  the  use  of  the  Gelfand-Levitan  algorithm.  To  be 
specific,  let  us  define  R(z)  by 

+°°  . , 

R(z)  =  (2tt) ” 1  /  b (k)  e"1KZ  dk.  (3.26) 

•.00 

The  reflection  coefficient  b(k)  satisfies  the  following  condi¬ 
tions  [see  Kay,  1954;  Kay  and  Moses,  1956;  Faddeev,  1958]. 

b(-k)  =  b*(k),  b(k)  analytic  in  the  upper  half-plane, 

b(0)  =  -1,  R(z)  e  0  for  z  <  2zQ.  (3.27) 

Let  us  define  the  Gelfand-Levitan  kernel  K(z,y)  by 
K(z,y)  =  0  for  either  y  >  z,  or  z  <  Zq.  (3.28) 

For  y  <  z  and  at  the  same  time  z  >  we  require  K(z,y)  to 
satisfy  the  Gelfand-Levitan  equation: 

z 

K(z,u)  R(u+y)  du 

(3.29) 

where  n(x)  is  the  Heaviside  function  defined  by  q (x)  =  1  for 
x  >  0,  n(x)  =  0  for  x  <  0. 

Having  found  the  Gelfand-Levitan  kernel  K(z,v),  V(z) 
is  given  by  the  simple  expression 

V (z)  =  2  K(z , z)  .  (3.30) 

However,  the  electric  field  E^(z)  can  also  ke  obtained  using 

E.  (z)  =  eikz  +  b(k)  e‘ikz  +  /  K(z,u)  [eiku  +  b(k)  e“iku]  du. 

k  2z0“z 

U  (3.31) 

Thus  to  find  V(z)  using  the  algorithm,  we  must  ob¬ 
tain  the  reflection  coefficient  b(k)  from  its  phase. 


K(z,y)  =  -  R(z+y)  -  n(z+y-2z«)  / 


This  is  accomplished  using  the  analytic  properties 
of  b(k)  and,  hence,  log  b(k)  and  a  generalized  form  of  the 
Hilbert  transform.  To  be  specific,  (d/dk)  <{>(k)  is  found  from 
the  time  delay  in  accordance  with  Eq.  (3.10).  The  phase  <f>(k) 
is  obtained  by  integration  with  the  boundary  condition  4>(0) 

=  7i .  Let  v(k)  be  defined  by 


v(k)  - 

4>  (k) 

-  2kzQ 

(3.32) 

and  w(k)  by 

w(k)  ■ 

log  1 

b  (k)  | 

(3.33) 

then 

w(k)  = 

-  T 

7T  J 

—  CO 

kv(k’) 

k1  (kT-k)  dk  • 

(3.34) 

In  Eq.  (3.34),  the  symbol  P  means  the  principal  part  of  the 
integral . 

It  should  be  mentioned  that  there  are  variational 
principles  available  [Moses,  19773  which  enable  one  to  obtain 
V(z)  when  b(k)  is  known.  These  principles  have  an  upper 
bound  built  into  them. 

One  also  has  available  a  generalization  of  the 
Celfand-Levitan  algorithm.  If,  for  a  given  reflection  coeffi¬ 
cient  bg(k),  one  knows  the  electron  density  Vq(z),  one  can  ob¬ 
tain  V(z)  -  Vq(z)  in  terms  of  b(k)  -  bg(k)  [Pechenick  and 
Cohen,  19813.  One  can  view  this  generalization  as  offering 
at  least  two  options.  One.  may  regard  V(z)  -  Vq(z)  as  the 
error  in  density  due  to  an  error  b(k)  -  bg(k)  in  the  reflec¬ 
tion  coefficient.  Or  one  may  think  of  Vq(z)  as  being  the  den¬ 
sity  associated  with  a  time  delay  leading  to  bg(k)  having  been 
obtained  from  a  model  or  a  previous  calculation  (even  using 
the  WKE  method).  Then  V(x)  is  obtained  as  a  relatively  small 
change  due  to  the  change  in  the  reflection  coefficient.  The 
variational  principle  can  also  be  used  to  obtain  V(z)  -  Vq(z) 
from  b(k)  -  bg(k)  together  with  a  bound  on  this  difference. 


3.2  Comparison  of  the  WKB  Method  with  the  Full-Wave  Method 


for  Profiles  for  Which  the  Full-Wave  Equation  can  be 

Solved  for  Exactly. 

3.2.1  General  Outline  of  Research  Procedure 

The  use  of  the  WKB  approximation  for  the  calculation 
of  electron  density  profiles  is  a  well-established  procedure. 
As  mentioned  in  the  previous  section,  the  WKB  approximation  is 
not  exact.  As  part  of  our  research  effort  we  are  considering 
potentials  V(z)  (which,  as  seen  above,  is  equivalent  to  con¬ 
sidering  number  densities  N(z)),  for  which  the  Schroedinger 
equation  can  be  solved  for  exactly  and  compare  the  exact 
time  delays  (or  equivalently  virtual  heights)  with  the  time 
delays  which  would  be  calculated  from  the  potential  using  the 
WKB  approximation  discussed  in  Section  3.1.  A  comparison  of 
the  time-delays  will  give  us  some  idea  of  how  good  the  WKB  ap¬ 
proximation  is.  Surprisingly,  perhaps,  our  work  in  this  area 
will  probably  be  one  of  the  few  to  compare  exact  time-delays 
with  WKB  time  delays,  even  though  the  WKB  time-delay  has  been 
used  for  almost  sixty  years.  The  reason  that  this  comparison 
was  not  made  earlier  is  that  the  Schroedinger  equation  had 
not  been  solved  exactly  for  potentials  V(x)  for  which  the  WFE 
approximation  might  be  expected  to  give  an  accurate  result. 
However,  in  more  recent  years  it  is  possible  to  construct  po¬ 
tentials  from  a  given  reflection  coefficient  for  which  the 
Schroedinger  equation  can  be  solved  exactly  and  in  terms  of 
elementary  functions.  Thus  having  specified  a  reflection  b(k) 
and  hence  as  part  of  the  reflection  coefficient,  its  phase  and 
thus  the  exact  time-delay,  we  compute  the  potential.  Having 
the  potential,  it  is  then  an  easy  matter  to  compute  the  WKB 
time-delay.  The  exact  time-delay  is  then  compared  with  the 
WKB  time -delay. 

Using  inverse  methods  it  is  possible  to  give  a 
class  of  potentials  which  are  smooth  and  continuous  and  which 
would  seem  to  model  a  single  layer  ionosphere  insofar  as  the 


general  shape  is  concerned.  Since  these  potentials  are  ob¬ 
tained  from  given  reflection  coefficients,  the  phase  and 
hence  the  time-delay  are  known  exactly.  From  such  examples 
one  can  compute  the  time-delay  using  the  UKB  method.  The  UKB 
time-delay  bears  no  relation  to  the  exact  time-delays.  In 
the  simplest  case,  which  we  discuss  below,  the  scaling  is  in¬ 
appropriate  for  the  UKE  method  to  be  valid  despite  the  ap¬ 
parent  smoothness.  Other  examples  remain  to  be  examined  for 
scaling . 


3.2.2  General  Properties  of  the  UKB  Time-Delay.  Failure 
of  the  UKB  Method.  Tentative  Conclusions. 

The  formula  for  the  UKB  time-delay  in  terms  of  a 
potential  is  identical  to  the  formula  for  obtaining  the  height 
of  a  hill  (essentially  V(x))  from  the  time  required  for  a  ball 
rolled  up  the  hill  with  a  given  velocity  to  come  down  again  to 
the  point  at  which  the  ball  started  rolling.  Generally,  the 
higher  the  initial  velocity  is,  the  greater  the  time  of  return 
will  be  because  the  ball  will  have  to  roll  higher  up  the  hill. 
In  wave  terms,  the  UKB  tine-delay  will  generally  increase  with 
increasing  k.  Our  examples  of  potentials  given  below  will 
show  that  the  UKB  behavior  for  the  time-delay  will  be  as  ex¬ 
pected.  However,  the  exact  time-delay  will  decrease  with  in¬ 
creasing  k  no  matter  how  slowly  growing  the  potential  of  the 
given  class  is.  In  the  simplest  case,  discussed  in  some  de¬ 
tail  below;  the  UKB  time-delay  bears  little  relationship  to 
the  exact,  full-wave  time-delay. 

Before  we  give  the  examples,  we  can  draw  some  conclu¬ 
sions  from  the  comparison  of  full-wave  and  UKB  theory.  It  is 
mathematically  possible  to  have  layers  for  which  the  UKB  time- 
delay  bears  no  relation  to  the  full-wave  time  delay.  In  prac¬ 
tice  one  usually  observes  time-delays  which  increase  with  k  in 
analogy  to  the  UKB  behavior.  Uhile  these  time-delays  may 
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arise  from  potentials  for  which  the  WKR  approximation  is  rea¬ 
sonably  good,  from  the  inverse  point  of  view  one  will  not  know 
that  the  observed  time-delay  (which  coincides  with  the  exact 
time  delay  for  the  potential)  can  be  used  in  a  WKB  approxima¬ 
tion  without  comparing  the  potential  from  the  I7KB  approxima¬ 
tion  with  that  of  the  full-wave  theory.  Since  it  is  necessary 
to  use  the  full  wave  theory  for  the  comparison,  and  since  in 
principle  and  probably  in  practice  (when  the  numerical  methods 
are  fully  established) ,  it  is  as  easy  to  use  the  full-wave  in¬ 
version  method  as  the  WKB  inversion,  one  should  use  the  exact 
method  in  the  first  place  to  prevent  error.  The  WK3  method 
may  also  be  used  afterward  if  one  chooses,  but  only  to  see 
how  well  it  predicts  the  exact  potential. 


3.3  Examples 

Our  examples  arise  through  the  choice  of  a  suitable 
set  of  reflection  coefficients  b(k).  Let  us  consider  the 
Gelfand-Levitan  equations  (3.26)  and  (3.29).  It  is  easy  to 
see  that  near  z  =  Zq,  V(z)  is  given  by 

V (z)  =  2  ^  (2tt)_1  /  b (k)  e"2ikz  dk.  (3.35) 


The  class  of  reflection  coefficients  we  shall  con¬ 
sider  are  dependent  on  two  parameters,  a  positive  constant  a, 
and  a  positive  integer  n: 


b  (k) 


(ia)n  # 
(k+ia)n 


(3.36) 


We  have  picked  b(k)  so  that  it  satisfies 

b (0)  =  -  1.  (3.37) 

This  condition  is  necessary  for  V(z)  to  be  positive.  All  the 
corresponding  potentials  Vn(z)  are  zero  for  negative  values 
of  z.  We  thus  are  taking  Zq  to  be  zero  and  regarding  the 
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transmitter  to  be  at  the  bottom  end  of  the  layer.  [The  gener¬ 
alization  to  positive  z q  simply  adds  a  time-delay  to  free- 
space  propagation  from  the  transmitter  at  z  =  0  to  the  bottom 
of  the  layer  at  z  =  Zq.]  It  is  also  seen  that  (in  the  vicin¬ 
ity  of  z  -  0)  V  (z)  is  of  the  form 

Vn(z)  =  C^z11-^  ,  z  greater  than  or  equal  to  zero,  (3.38) 

for  n  greater  than  or  equal  to  2  where  is  a  real  positive 
constant.  For  n  =  2  we  see  from  (3.35)  that  the  potential 
has  a  jump  at  z  =  0  from  zero  for  z  less  than  zero  to  C2  for  z 
greater  than  but  near  zero.  For  n  =  3  the  potential  is  con¬ 
tinuous  at  zero  and  ^3(0)  =  0.  For  n  =  A,  V  and  the  first 
derivative  of  Vn  are  zero  at  z  =  0  and  are  continuous  there. 

As  one  takes  higher  values  of  n  an  increasing  number  of  deriv¬ 
atives  are  continuous  and  zero  at  z  =  0.  Thus  one  can  make 
the  potential  increase  from  zero  to  finite  values  as  slowly  as 
one  wishes.  High  values  of  n  correspond  to  increasingly 
smoother  variations  of  V  from  zero  to  positive  values.  In 
regions  sufficiently  near  z  =  0,  where  varies  in  as  smooth 
a  manner  as  one  wishes,  one  might  expect  the  best  chances  of 
the  \TKB  time-delay  to  be  valid.  However,  for  all  values  of  k 
and  n  the  exact  time  delay  decreases ,  as  we  shall  show  shortly 

A  possible  explanation  for  the  failure  of  the  UF.B 
method  is  the  fact  that  in  addition  to  requiring  that  the  in¬ 
dex  of  refraction  vary  slowly  and  not  be  too  small,  the  value 
of  k  should  not  be  small  [Budden,  1966],  For  low  value  of  n, 
k  is  small  compared  to  the  scale  a,  but  for  high  values  of  n, 
k  can  be  made  increasingly  large  and  the  blCB  procedure  goes 
through  formally  but  still  gives  the  wrong  character  for  the 
tine-delay  as  compared  with  the  exact  time-delay.  These  re¬ 
sults  suggest  more  strongly  than  ever  that  the  WKF  inverse 
method  requires  a  closer  scrutiny  even  where  it  is  believed 
to  be  valid. 
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The  time-delay  Tn(k)  arising  from  the  reflection  co¬ 
efficient  of  Eq.  (3.36)  is  readily  found  to  be 


TOO  = 

n  2 


na 


k2+a‘ 


(3.39) 


The  exact  potential  Vn(z)  have  been  computed  and 
sketched  by  Pechenick  and  Cohen  [1981].  As  n  increases  the 
potentials,  which  may  be  considered  to  be  a  two-parameter  fam¬ 
ily  of  functions  depending  on  n  and  a,  are  all  positive  ar.d, 
for  n  equal  to  or  greater  than  3  start  at  zero  for  z  =  0  andi 
have  a  suitable  number  of  derivatives  vanishing  at  z  =  0. 

The  potentials  increase  monotonically ,  reach  a  maximum  and 
then  decay  exponentially  for  large  z.  For  a  given  value  of  a, 
higher  values  of  n  correspond  to  flatter  potentials,  i.e.  the 
effective  base  of  the  potential  increases  with  respect  to  the 
height . 

In  Figure  19,  we  have  plotted  To  the  eye  it 

looks  like  a  "normal"  potential.  One  would  think  that  the 
VJKB  calculation  of  the  time  delay  for  that  portion  of  the  po¬ 
tential  for  z  between  0  and  the  value  of  z  corresponding  to 
the  maximum  would  be  a  good  approximation.  In  Figure  20,  the 
exact  and  WKB  values  of  the  time-delay  are  plotted  as  a  func¬ 
tion  of  k.  The  discrepancy  between  the  two  time  delays  is 
astonishing!  The  WKB  time-delay  does  not  even  approach  the 
values  of  the  exact  one  for  most  values  of  k.  Thus,  even  for 
potentials  which  look  as  though  VJKB  time  delay  ought  to  be 
valid,  this  is  not  generally  the  case. 

It  should  be  noted  that  V^(x)  does  not  model  the 
ionosphere  at  all  well  because  of  the  length  to  height  ratio, 
which  indeed  precludes  the  use  of  the  WKB  method.  A  better 
exact  model,  which  we  are  trying  to  obtain,  will  be  more  use¬ 
ful  for  comparing  the  WKB  and  exact  methods. 


Figure  20.  Full  Wave  and  WKB  Virtual  Heights 


3.4  Future  Research  Directions 

Having  seen  that  potentials  exist  to  which  the  WKB 
approximation  can  not  be  applied,  we  intend  to  introduce  po¬ 
tentials  for  which  the  WKB  time-delay  has  some  of  the  proper¬ 
ties  of  the  exact  time-delay  to  find  the  domain  where  WKB  has 
meaning. 

In  addition  to  obtaining  potentials  for  which  the 
WKB  time-delay  makes  some  sense,  we  shall  give  examples  for 
even  odder  behavior  of  the  exact  time-delay  than  those  of  the 
present  report.  They  will  include  potentials  for  which  the 
time-delay  is  zero  for  all  k  and  potentials  for  which  the 
time-delay  is  negative:  causality  seems  to  be  violated.  The 
potentials  for  these  unexpected  time-delays  look  "normal". 

The  result  of  these  investigations  appear  to  indi¬ 
cate  that  in  many  circumstances  time-delay  is  not  a  particu¬ 
larly  useful  notion.  More  useful  will  be  the  shape  of  the 
returned  wave  for  which  the  shape  of  the  pulse  is  known.  By 
deconvolution  techniques  the  reflection  coefficient  b(k)  can 
be  found  and  one  can  then  use  the  full-wave  solutions.  These 
possibilities  will  also  be  examined  from  a  numerical  point  of 
view. 
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